/* -*- c -*- */

/* The purpose of this module is to add faster math for array scalars
   that does not go through the ufunc machinery

   but still supports error-modes.
*/
#define PY_SSIZE_T_CLEAN
#include <Python.h>

#define _UMATHMODULE
#define _MULTIARRAYMODULE
#define NPY_NO_DEPRECATED_API NPY_API_VERSION

#include "npy_config.h"
#include "numpy/arrayobject.h"
#include "numpy/ufuncobject.h"
#include "numpy/arrayscalars.h"

#include "npy_import.h"
#include "npy_pycompat.h"

#include "numpy/halffloat.h"
#include "templ_common.h"

#include "binop_override.h"
#include "npy_longdouble.h"

/* Basic operations:
 *
 *  BINARY:
 *
 * add, subtract, multiply, divide, remainder, divmod, power,
 * floor_divide, true_divide
 *
 * lshift, rshift, and, or, xor (integers only)
 *
 * UNARY:
 *
 * negative, positive, absolute, nonzero, invert, int, long, float, oct, hex
 *
 */

/**begin repeat
 *  #name = byte, short, int, long, longlong#
 *  #type = npy_byte, npy_short, npy_int, npy_long, npy_longlong#
 */
static void
@name@_ctype_add(@type@ a, @type@ b, @type@ *out) {
    *out = a + b;
    if ((*out^a) >= 0 || (*out^b) >= 0) {
        return;
    }
    npy_set_floatstatus_overflow();
    return;
}
static void
@name@_ctype_subtract(@type@ a, @type@ b, @type@ *out) {
    *out = a - b;
    if ((*out^a) >= 0 || (*out^~b) >= 0) {
        return;
    }
    npy_set_floatstatus_overflow();
    return;
}
/**end repeat**/

/**begin repeat
 *  #name = ubyte, ushort, uint, ulong, ulonglong#
 *  #type = npy_ubyte, npy_ushort, npy_uint, npy_ulong, npy_ulonglong#
 */
static void
@name@_ctype_add(@type@ a, @type@ b, @type@ *out) {
    *out = a + b;
    if (*out >= a && *out >= b) {
        return;
    }
    npy_set_floatstatus_overflow();
    return;
}
static void
@name@_ctype_subtract(@type@ a, @type@ b, @type@ *out) {
    *out = a - b;
    if (a >= b) {
        return;
    }
    npy_set_floatstatus_overflow();
    return;
}
/**end repeat**/

#ifndef NPY_SIZEOF_BYTE
#define NPY_SIZEOF_BYTE 1
#endif

/**begin repeat
 *
 * #name = byte, ubyte, short, ushort,
 *         int, uint, long, ulong#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort,
 *         npy_int, npy_uint, npy_long, npy_ulong#
 * #big =  npy_int, npy_uint, npy_int, npy_uint,
 *         npy_longlong, npy_ulonglong, npy_longlong, npy_ulonglong#
 * #NAME = BYTE, UBYTE, SHORT, USHORT,
 *         INT, UINT, LONG, ULONG#
 * #SIZENAME = BYTE*2, SHORT*2, INT*2, LONG*2#
 * #SIZE = INT*4,LONGLONG*4#
 * #neg = (1,0)*4#
 */
#if NPY_SIZEOF_@SIZE@ > NPY_SIZEOF_@SIZENAME@
static void
@name@_ctype_multiply(@type@ a, @type@ b, @type@ *out) {
    @big@ temp;
    temp = ((@big@) a) * ((@big@) b);
    *out = (@type@) temp;
#if @neg@
    if (temp > NPY_MAX_@NAME@ || temp < NPY_MIN_@NAME@)
#else
        if (temp > NPY_MAX_@NAME@)
#endif
            npy_set_floatstatus_overflow();
    return;
}
#endif
/**end repeat**/

/**begin repeat
 *
 * #name = int, uint, long, ulong,
 *         longlong, ulonglong#
 * #type = npy_int, npy_uint, npy_long, npy_ulong,
 *         npy_longlong, npy_ulonglong#
 * #SIZE = INT*2, LONG*2, LONGLONG*2#
 */
#if NPY_SIZEOF_LONGLONG == NPY_SIZEOF_@SIZE@
static void
@name@_ctype_multiply(@type@ a, @type@ b, @type@ *out) {
    if (npy_mul_with_overflow_@name@(out, a, b)) {
        npy_set_floatstatus_overflow();
    }
    return;
}
#endif
/**end repeat**/

/**begin repeat
 *
 * #name = byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong#
 * #neg = (1,0)*5#
 */
static void
@name@_ctype_divide(@type@ a, @type@ b, @type@ *out) {
    if (b == 0) {
        npy_set_floatstatus_divbyzero();
        *out = 0;
    }
#if @neg@
    else if (b == -1 && a < 0 && a == -a) {
        npy_set_floatstatus_overflow();
        *out = a / b;
    }
#endif
    else {
#if @neg@
        @type@ tmp;
        tmp = a / b;
        if (((a > 0) != (b > 0)) && (a % b != 0)) {
            tmp--;
        }
        *out = tmp;
#else
        *out = a / b;
#endif
    }
}

#define @name@_ctype_floor_divide @name@_ctype_divide

static void
@name@_ctype_remainder(@type@ a, @type@ b, @type@ *out) {
    if (a == 0 || b == 0) {
        if (b == 0) npy_set_floatstatus_divbyzero();
        *out = 0;
        return;
    }
#if @neg@
    else if ((a > 0) == (b > 0)) {
        *out = a % b;
    }
    else {
        /* handled like Python does */
        *out = a % b;
        if (*out) *out += b;
    }
#else
    *out = a % b;
#endif
}
/**end repeat**/

/**begin repeat
 *
 * #name = byte, ubyte, short, ushort, int, uint, long,
 *         ulong, longlong, ulonglong#
 * #otyp = npy_float*4, npy_double*6#
 */
#define @name@_ctype_true_divide(a, b, out)     \
    *(out) = ((@otyp@) (a)) / ((@otyp@) (b));
/**end repeat**/

/* b will always be positive in this call */
/**begin repeat
 *
 * #name = byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong#
 * #upc = BYTE, UBYTE, SHORT, USHORT, INT, UINT,
 *        LONG, ULONG, LONGLONG, ULONGLONG#
 */
static void
@name@_ctype_power(@type@ a, @type@ b, @type@ *out) {
    @type@ tmp;

    if (b == 0) {
        *out = 1;
        return;
    }
    if (a == 1) {
        *out = 1;
        return;
    }

    tmp = b & 1 ? a : 1;
    b >>= 1;
    while (b > 0) {
        a *= a;
        if (b & 1) {
            tmp *= a;
        }
        b >>= 1;
    }
    *out = tmp;
}
/**end repeat**/


/**begin repeat
 * #name = byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong#
 * #suffix = hh,uhh,h,uh,,u,l,ul,ll,ull#
 */

/**begin repeat1
 * #oper = and, xor, or#
 * #op = &, ^, |#
 */

#define @name@_ctype_@oper@(arg1, arg2, out) *(out) = (arg1) @op@ (arg2)

/**end repeat1**/

#define @name@_ctype_lshift(arg1, arg2, out) *(out) = npy_lshift@suffix@(arg1, arg2)
#define @name@_ctype_rshift(arg1, arg2, out) *(out) = npy_rshift@suffix@(arg1, arg2)

/**end repeat**/

/**begin repeat
 * #name = float, double, longdouble#
 * #type = npy_float, npy_double, npy_longdouble#
 * #c = f, , l#
 */
#define @name@_ctype_add(a, b, outp) *(outp) = (a) + (b)
#define @name@_ctype_subtract(a, b, outp) *(outp) = (a) - (b)
#define @name@_ctype_multiply(a, b, outp) *(outp) = (a) * (b)
#define @name@_ctype_divide(a, b, outp) *(outp) = (a) / (b)
#define @name@_ctype_true_divide @name@_ctype_divide


static void
@name@_ctype_floor_divide(@type@ a, @type@ b, @type@ *out) {
    *out = npy_floor_divide@c@(a, b);
}


static void
@name@_ctype_remainder(@type@ a, @type@ b, @type@ *out) {
    *out = npy_remainder@c@(a, b);
}


static void
@name@_ctype_divmod(@type@ a, @type@ b, @type@ *out1, @type@ *out2) {
    *out1 = npy_divmod@c@(a, b, out2);
}


/**end repeat**/

#define half_ctype_add(a, b, outp) *(outp) = \
        npy_float_to_half(npy_half_to_float(a) + npy_half_to_float(b))
#define half_ctype_subtract(a, b, outp) *(outp) = \
        npy_float_to_half(npy_half_to_float(a) - npy_half_to_float(b))
#define half_ctype_multiply(a, b, outp) *(outp) = \
        npy_float_to_half(npy_half_to_float(a) * npy_half_to_float(b))
#define half_ctype_divide(a, b, outp) *(outp) = \
        npy_float_to_half(npy_half_to_float(a) / npy_half_to_float(b))
#define half_ctype_true_divide half_ctype_divide


static void
half_ctype_floor_divide(npy_half a, npy_half b, npy_half *out) {
    npy_half mod;

    if (!b) {
        *out = a / b;
    } else {
        *out = npy_half_divmod(a, b, &mod);
    }
}


static void
half_ctype_remainder(npy_half a, npy_half b, npy_half *out) {
    npy_half_divmod(a, b, out);
}


static void
half_ctype_divmod(npy_half a, npy_half b, npy_half *out1, npy_half *out2) {
    *out1 = npy_half_divmod(a, b, out2);
}

/**begin repeat
 * #name = cfloat, cdouble, clongdouble#
 * #rname = float, double, longdouble#
 * #rtype = npy_float, npy_double, npy_longdouble#
 * #c = f,,l#
 */
#define @name@_ctype_add(a, b, outp) do{        \
    (outp)->real = (a).real + (b).real;         \
    (outp)->imag = (a).imag + (b).imag;         \
    } while(0)
#define @name@_ctype_subtract(a, b, outp) do{   \
    (outp)->real = (a).real - (b).real;         \
    (outp)->imag = (a).imag - (b).imag;         \
    } while(0)
#define @name@_ctype_multiply(a, b, outp) do{                   \
    (outp)->real = (a).real * (b).real - (a).imag * (b).imag;   \
    (outp)->imag = (a).real * (b).imag + (a).imag * (b).real;   \
    } while(0)
/* Algorithm identical to that in loops.c.src, for consistency */
#define @name@_ctype_divide(a, b, outp) do{                         \
    @rtype@ in1r = (a).real;                                        \
    @rtype@ in1i = (a).imag;                                        \
    @rtype@ in2r = (b).real;                                        \
    @rtype@ in2i = (b).imag;                                        \
    @rtype@ in2r_abs = npy_fabs@c@(in2r);                           \
    @rtype@ in2i_abs = npy_fabs@c@(in2i);                           \
    if (in2r_abs >= in2i_abs) {                                     \
        if (in2r_abs == 0 && in2i_abs == 0) {                       \
            /* divide by zero should yield a complex inf or nan */  \
            (outp)->real = in1r/in2r_abs;                           \
            (outp)->imag = in1i/in2i_abs;                           \
        }                                                           \
        else {                                                      \
            @rtype@ rat = in2i/in2r;                                \
            @rtype@ scl = 1.0@c@/(in2r + in2i*rat);                 \
            (outp)->real = (in1r + in1i*rat)*scl;                   \
            (outp)->imag = (in1i - in1r*rat)*scl;                   \
        }                                                           \
    }                                                               \
    else {                                                          \
        @rtype@ rat = in2r/in2i;                                    \
        @rtype@ scl = 1.0@c@/(in2i + in2r*rat);                     \
        (outp)->real = (in1r*rat + in1i)*scl;                       \
        (outp)->imag = (in1i*rat - in1r)*scl;                       \
    }                                                               \
    } while(0)

#define @name@_ctype_true_divide @name@_ctype_divide

#define @name@_ctype_floor_divide(a, b, outp) do {      \
    @rname@_ctype_floor_divide(                         \
        ((a).real*(b).real + (a).imag*(b).imag),        \
        ((b).real*(b).real + (b).imag*(b).imag),        \
        &((outp)->real));                               \
    (outp)->imag = 0;                                   \
    } while(0)
/**end repeat**/



/**begin repeat
 * #name = byte, ubyte, short, ushort, int, uint, long, ulong,
 *         longlong, ulonglong, cfloat, cdouble, clongdouble#
 */
#define @name@_ctype_divmod(a, b, out, out2) {  \
    @name@_ctype_floor_divide(a, b, out);       \
    @name@_ctype_remainder(a, b, out2);         \
    }
/**end repeat**/


/**begin repeat
 * #name = float, double, longdouble#
 * #type = npy_float, npy_double, npy_longdouble#
 * #c = f,,l#
 */

static void
@name@_ctype_power(@type@ a, @type@ b, @type@ *out)
{
    *out = npy_pow@c@(a, b);
}

/**end repeat**/
static void
half_ctype_power(npy_half a, npy_half b, npy_half *out)
{
    const npy_float af = npy_half_to_float(a);
    const npy_float bf = npy_half_to_float(b);
    const npy_float outf = npy_powf(af,bf);
    *out = npy_float_to_half(outf);
}

/**begin repeat
 * #name = byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong,
 *         float, double, longdouble#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *         npy_float, npy_double, npy_longdouble#
 * #uns = (0,1)*5,0*3#
 */
static void
@name@_ctype_negative(@type@ a, @type@ *out)
{
#if @uns@
    npy_set_floatstatus_overflow();
#endif
    *out = -a;
}
/**end repeat**/

static void
half_ctype_negative(npy_half a, npy_half *out)
{
    *out = a^0x8000u;
}


/**begin repeat
 * #name = cfloat, cdouble, clongdouble#
 * #type = npy_cfloat, npy_cdouble, npy_clongdouble#
 */
static void
@name@_ctype_negative(@type@ a, @type@ *out)
{
    out->real = -a.real;
    out->imag = -a.imag;
}
/**end repeat**/

/**begin repeat
 * #name = byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong,
 *         half, float, double, longdouble#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *         npy_half, npy_float, npy_double, npy_longdouble#
 */
static void
@name@_ctype_positive(@type@ a, @type@ *out)
{
    *out = a;
}
/**end repeat**/

/**begin repeat
 * #name = cfloat, cdouble, clongdouble#
 * #type = npy_cfloat, npy_cdouble, npy_clongdouble#
 * #c = f,,l#
 */
static void
@name@_ctype_positive(@type@ a, @type@ *out)
{
    out->real = a.real;
    out->imag = a.imag;
}

static void
@name@_ctype_power(@type@ a, @type@ b, @type@ *out)
{
    *out = npy_cpow@c@(a, b);
}
/**end repeat**/


/**begin repeat
 * #name = ubyte, ushort, uint, ulong, ulonglong#
 */

#define @name@_ctype_absolute @name@_ctype_positive

/**end repeat**/


/**begin repeat
 * #name = byte, short, int, long, longlong#
 * #type = npy_byte, npy_short, npy_int, npy_long, npy_longlong#
 */
static void
@name@_ctype_absolute(@type@ a, @type@ *out)
{
    *out = (a < 0 ? -a : a);
}
/**end repeat**/

/**begin repeat
 * #name = float, double, longdouble#
 * #type = npy_float, npy_double, npy_longdouble#
 * #c = f,,l#
 */
static void
@name@_ctype_absolute(@type@ a, @type@ *out)
{
    *out = npy_fabs@c@(a);
}
/**end repeat**/

static void
half_ctype_absolute(npy_half a, npy_half *out)
{
    *out = a&0x7fffu;
}

/**begin repeat
 * #name = cfloat, cdouble, clongdouble#
 * #type = npy_cfloat, npy_cdouble, npy_clongdouble#
 * #rtype = npy_float, npy_double, npy_longdouble#
 * #c = f,,l#
 */
static void
@name@_ctype_absolute(@type@ a, @rtype@ *out)
{
    *out = npy_cabs@c@(a);
}
/**end repeat**/

/**begin repeat
 * #name = byte, ubyte, short, ushort, int, uint, long,
 *         ulong, longlong, ulonglong#
 */

#define @name@_ctype_invert(a, out) *(out) = ~a;

/**end repeat**/

/*** END OF BASIC CODE **/


/* The general strategy for commutative binary operators is to
 *
 * 1) Convert the types to the common type if both are scalars (0 return)
 * 2) If both are not scalars use ufunc machinery (-2 return)
 * 3) If both are scalars but cannot be cast to the right type
 * return NotImplemented (-1 return)
 *
 * 4) Perform the function on the C-type.
 * 5) If an error condition occurred, check to see
 * what the current error-handling is and handle the error.
 *
 * 6) Construct and return the output scalar.
 */

/**begin repeat
 * #name = byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong,
 *         half, float, longdouble,
 *         cfloat, cdouble, clongdouble#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *         npy_half, npy_float, npy_longdouble,
 *         npy_cfloat, npy_cdouble, npy_clongdouble#
 * #Name = Byte, UByte, Short, UShort, Int, UInt,
 *         Long, ULong, LongLong, ULongLong,
 *         Half, Float, LongDouble,
 *         CFloat, CDouble, CLongDouble#
 * #TYPE = NPY_BYTE, NPY_UBYTE, NPY_SHORT, NPY_USHORT, NPY_INT, NPY_UINT,
 *         NPY_LONG, NPY_ULONG, NPY_LONGLONG, NPY_ULONGLONG,
 *         NPY_HALF, NPY_FLOAT, NPY_LONGDOUBLE,
 *         NPY_CFLOAT, NPY_CDOUBLE, NPY_CLONGDOUBLE#
 */

static int
_@name@_convert_to_ctype(PyObject *a, @type@ *arg1)
{
    PyObject *temp;

    if (PyArray_IsScalar(a, @Name@)) {
        *arg1 = PyArrayScalar_VAL(a, @Name@);
        return 0;
    }
    else if (PyArray_IsScalar(a, Generic)) {
        PyArray_Descr *descr1;

        if (!PyArray_IsScalar(a, Number)) {
            return -1;
        }
        descr1 = PyArray_DescrFromTypeObject((PyObject *)Py_TYPE(a));
        if (PyArray_CanCastSafely(descr1->type_num, @TYPE@)) {
            PyArray_CastScalarDirect(a, descr1, arg1, @TYPE@);
            Py_DECREF(descr1);
            return 0;
        }
        else {
            Py_DECREF(descr1);
            return -1;
        }
    }
    else if (PyArray_GetPriority(a, NPY_PRIORITY) > NPY_PRIORITY) {
        return -2;
    }
    else if ((temp = PyArray_ScalarFromObject(a)) != NULL) {
        int retval = _@name@_convert_to_ctype(temp, arg1);

        Py_DECREF(temp);
        return retval;
    }
    return -2;
}

/**end repeat**/


/* Same as above but added exact checks against known python types for speed */

/**begin repeat
 * #name = double#
 * #type = npy_double#
 * #Name = Double#
 * #TYPE = NPY_DOUBLE#
 * #PYCHECKEXACT = PyFloat_CheckExact#
 * #PYEXTRACTCTYPE = PyFloat_AS_DOUBLE#
 */

static int
_@name@_convert_to_ctype(PyObject *a, @type@ *arg1)
{
    PyObject *temp;

    if (@PYCHECKEXACT@(a)){
        *arg1 = @PYEXTRACTCTYPE@(a);
        return 0;
    }

    if (PyArray_IsScalar(a, @Name@)) {
        *arg1 = PyArrayScalar_VAL(a, @Name@);
        return 0;
    }
    else if (PyArray_IsScalar(a, Generic)) {
        PyArray_Descr *descr1;

        if (!PyArray_IsScalar(a, Number)) {
            return -1;
        }
        descr1 = PyArray_DescrFromTypeObject((PyObject *)Py_TYPE(a));
        if (PyArray_CanCastSafely(descr1->type_num, @TYPE@)) {
            PyArray_CastScalarDirect(a, descr1, arg1, @TYPE@);
            Py_DECREF(descr1);
            return 0;
        }
        else {
            Py_DECREF(descr1);
            return -1;
        }
    }
    else if (PyArray_GetPriority(a, NPY_PRIORITY) > NPY_PRIORITY) {
        return -2;
    }
    else if ((temp = PyArray_ScalarFromObject(a)) != NULL) {
        int retval = _@name@_convert_to_ctype(temp, arg1);

        Py_DECREF(temp);
        return retval;
    }
    return -2;
}

/**end repeat**/


/**begin repeat
 * #name = byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong,
 *         half, float, double, cfloat, cdouble#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *         npy_half, npy_float, npy_double, npy_cfloat, npy_cdouble#
 */
static int
_@name@_convert2_to_ctypes(PyObject *a, @type@ *arg1,
                           PyObject *b, @type@ *arg2)
{
    int ret;
    ret = _@name@_convert_to_ctype(a, arg1);
    if (ret < 0) {
        return ret;
    }
    ret = _@name@_convert_to_ctype(b, arg2);
    if (ret < 0) {
        return ret;
    }
    return 0;
}
/**end repeat**/

/**begin repeat
 * #name = longdouble, clongdouble#
 * #type = npy_longdouble, npy_clongdouble#
 */

static int
_@name@_convert2_to_ctypes(PyObject *a, @type@ *arg1,
                           PyObject *b, @type@ *arg2)
{
    int ret;
    ret = _@name@_convert_to_ctype(a, arg1);
    if (ret == -2) {
        ret = -3;
    }
    if (ret < 0) {
        return ret;
    }
    ret = _@name@_convert_to_ctype(b, arg2);
    if (ret == -2) {
        ret = -3;
    }
    if (ret < 0) {
        return ret;
    }
    return 0;
}

/**end repeat**/


/**begin repeat
 *
 * #name = (byte, ubyte, short, ushort, int, uint,
 *             long, ulong, longlong, ulonglong)*12,
 *         (half, float, double, longdouble,
 *             cfloat, cdouble, clongdouble)*5,
 *         (half, float, double, longdouble)*2#
 * #Name = (Byte, UByte, Short, UShort, Int, UInt,
 *             Long, ULong,LongLong,ULongLong)*12,
 *         (Half, Float, Double, LongDouble,
 *             CFloat, CDouble, CLongDouble)*5,
 *         (Half, Float, Double, LongDouble)*2#
 * #type = (npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *             npy_long, npy_ulong, npy_longlong, npy_ulonglong)*12,
 *         (npy_half, npy_float, npy_double, npy_longdouble,
 *             npy_cfloat, npy_cdouble, npy_clongdouble)*5,
 *         (npy_half, npy_float, npy_double, npy_longdouble)*2#
 *
 * #oper = add*10, subtract*10, multiply*10, remainder*10,
 *         divmod*10, floor_divide*10, lshift*10, rshift*10, and*10,
 *         or*10, xor*10, true_divide*10,
 *         add*7, subtract*7, multiply*7, floor_divide*7, true_divide*7,
 *         divmod*4, remainder*4#
 *
 * #fperr = 1*60,0*50,1*10,
 *          1*35,
 *          1*8#
 * #twoout = 0*40,1*10,0*70,
 *           0*35,
 *           1*4,0*4#
 * #otype = (npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *             npy_long, npy_ulong, npy_longlong, npy_ulonglong)*11,
 *         npy_float*4, npy_double*6,
 *         (npy_half, npy_float, npy_double, npy_longdouble,
 *             npy_cfloat, npy_cdouble, npy_clongdouble)*5,
 *         (npy_half, npy_float, npy_double, npy_longdouble)*2#
 * #OName = (Byte, UByte, Short, UShort, Int, UInt,
 *              Long, ULong, LongLong, ULongLong)*11,
 *          Float*4, Double*6,
 *          (Half, Float, Double, LongDouble,
 *              CFloat, CDouble, CLongDouble)*5,
 *          (Half, Float, Double, LongDouble)*2#
 */

static PyObject *
@name@_@oper@(PyObject *a, PyObject *b)
{
    PyObject *ret;
    @type@ arg1, arg2;
    @otype@ out;

#if @twoout@
    @otype@ out2;
    PyObject *obj;
#endif

#if @fperr@
    int retstatus;
    int first;
#endif

    BINOP_GIVE_UP_IF_NEEDED(a, b, nb_@oper@, @name@_@oper@);

    switch(_@name@_convert2_to_ctypes(a, &arg1, b, &arg2)) {
        case 0:
            break;
        case -1:
            /* one of them can't be cast safely must be mixed-types*/
            return PyArray_Type.tp_as_number->nb_@oper@(a,b);
        case -2:
            /* use default handling */
            if (PyErr_Occurred()) {
                return NULL;
            }
            return PyGenericArrType_Type.tp_as_number->nb_@oper@(a,b);
        case -3:
            /*
             * special case for longdouble and clongdouble
             * because they have a recursive getitem in their dtype
             */
            Py_INCREF(Py_NotImplemented);
            return Py_NotImplemented;
    }

#if @fperr@
    npy_clear_floatstatus_barrier((char*)&out);
#endif

    /*
     * here we do the actual calculation with arg1 and arg2
     * as a function call.
     */
#if @twoout@
    @name@_ctype_@oper@(arg1, arg2, (@otype@ *)&out, &out2);
#else
    @name@_ctype_@oper@(arg1, arg2, (@otype@ *)&out);
#endif

#if @fperr@
    /* Check status flag.  If it is set, then look up what to do */
    retstatus = npy_get_floatstatus_barrier((char*)&out);
    if (retstatus) {
        int bufsize, errmask;
        PyObject *errobj;

        if (PyUFunc_GetPyValues("@name@_scalars", &bufsize, &errmask,
                                &errobj) < 0) {
            return NULL;
        }
        first = 1;
        if (PyUFunc_handlefperr(errmask, errobj, retstatus, &first)) {
            Py_XDECREF(errobj);
            return NULL;
        }
        Py_XDECREF(errobj);
    }
#endif


#if @twoout@
    ret = PyTuple_New(2);
    if (ret == NULL) {
        return NULL;
    }
    obj = PyArrayScalar_New(@OName@);
    if (obj == NULL) {
        Py_DECREF(ret);
        return NULL;
    }
    PyArrayScalar_ASSIGN(obj, @OName@, out);
    PyTuple_SET_ITEM(ret, 0, obj);
    obj = PyArrayScalar_New(@OName@);
    if (obj == NULL) {
        Py_DECREF(ret);
        return NULL;
    }
    PyArrayScalar_ASSIGN(obj, @OName@, out2);
    PyTuple_SET_ITEM(ret, 1, obj);
#else
    ret = PyArrayScalar_New(@OName@);
    if (ret == NULL) {
        return NULL;
    }
    PyArrayScalar_ASSIGN(ret, @OName@, out);
#endif
    return ret;
}

/**end repeat**/

#define _IS_ZERO(x) (x == 0)

/**begin repeat
 *
 * #name = byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong,
 *         half, float, double, longdouble,
 *         cfloat, cdouble, clongdouble#
 *
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *         npy_half, npy_float, npy_double, npy_longdouble,
 *         npy_cfloat, npy_cdouble, npy_clongdouble#
 *
 * #Name = Byte, UByte, Short, UShort, Int, UInt,
 *         Long, ULong, LongLong, ULongLong,
 *         Half, Float, Double, LongDouble,
 *         CFloat, CDouble, CLongDouble#
 *
 * #oname = float*4, double*6, half, float, double, longdouble,
 *          cfloat, cdouble, clongdouble#
 *
 * #otype = npy_float*4, npy_double*6, npy_half, npy_float,
 *          npy_double, npy_longdouble,
 *          npy_cfloat, npy_cdouble, npy_clongdouble#
 *
 * #OName = Float*4, Double*6, Half, Float,
 *          Double, LongDouble,
 *          CFloat, CDouble, CLongDouble#
 *
 * #isint = 1*10,0*7#
 * #isuint = (0,1)*5,0*7#
 * #cmplx = 0*14,1*3#
 * #iszero = _IS_ZERO*10, npy_half_iszero, _IS_ZERO*6#
 * #zero = 0*10, NPY_HALF_ZERO, 0*6#
 * #one = 1*10, NPY_HALF_ONE, 1*6#
 */

static PyObject *
@name@_power(PyObject *a, PyObject *b, PyObject *modulo)
{
    PyObject *ret;
    @type@ arg1, arg2, out;

    BINOP_GIVE_UP_IF_NEEDED(a, b, nb_power, @name@_power);

    switch(_@name@_convert2_to_ctypes(a, &arg1, b, &arg2)) {
        case 0:
            break;
        case -1:
            /* can't cast both safely mixed-types? */
            return PyArray_Type.tp_as_number->nb_power(a,b,modulo);
        case -2:
            /* use default handling */
            if (PyErr_Occurred()) {
                return NULL;
            }
            return PyGenericArrType_Type.tp_as_number->nb_power(a,b,modulo);
        case -3:
        default:
            /*
             * special case for longdouble and clongdouble
             * because they have a recursive getitem in their dtype
             */
            Py_INCREF(Py_NotImplemented);
            return Py_NotImplemented;
    }

    if (modulo != Py_None) {
        /* modular exponentiation is not implemented (gh-8804) */
        Py_INCREF(Py_NotImplemented);
        return Py_NotImplemented;
    }

#if !@isint@
    npy_clear_floatstatus_barrier((char*)&out);
#endif
    /*
     * here we do the actual calculation with arg1 and arg2
     * as a function call.
     */
#if @isint@ && !@isuint@
    if (arg2 < 0) {
        PyErr_SetString(PyExc_ValueError,
                "Integers to negative integer powers are not allowed.");
        return NULL;
    }
#endif
    @name@_ctype_power(arg1, arg2, &out);

#if !@isint@
    /* Check status flag.  If it is set, then look up what to do */
    int retstatus = npy_get_floatstatus_barrier((char*)&out);
    if (retstatus) {
        int bufsize, errmask;
        PyObject *errobj;

        if (PyUFunc_GetPyValues("@name@_scalars", &bufsize, &errmask,
                                &errobj) < 0) {
            return NULL;
        }
        int first = 1;
        if (PyUFunc_handlefperr(errmask, errobj, retstatus, &first)) {
            Py_XDECREF(errobj);
            return NULL;
        }
        Py_XDECREF(errobj);
    }
#endif

    ret = PyArrayScalar_New(@Name@);
    if (ret == NULL) {
        return NULL;
    }
    PyArrayScalar_ASSIGN(ret, @Name@, out);

    return ret;
}


/**end repeat**/
#undef _IS_ZERO


/**begin repeat
 *
 * #name = cfloat, cdouble#
 *
 */

/**begin repeat1
 *
 * #oper = divmod, remainder#
 *
 */

#define @name@_@oper@ NULL

/**end repeat1**/

/**end repeat**/

/**begin repeat
 *
 * #oper = divmod, remainder#
 *
 */

/* 
Complex numbers do not support remainder operations. Unfortunately, 
the type inference for long doubles is complicated, and if a remainder 
operation is not defined - if the relevant field is left NULL - then 
operations between long doubles and objects lead to an infinite recursion 
instead of a TypeError. This should ensure that once everything gets
converted to complex long doubles you correctly get a reasonably
informative TypeError. This fixes the last part of bug gh-18548.
*/

static PyObject *
clongdouble_@oper@(PyObject *a, PyObject *b)
{
    npy_clongdouble arg1, arg2;

    BINOP_GIVE_UP_IF_NEEDED(a, b, nb_@oper@, clongdouble_@oper@);

    switch(_clongdouble_convert2_to_ctypes(a, &arg1, b, &arg2)) {
        case 0:
            break;
        case -1:
            /* one of them can't be cast safely must be mixed-types*/
            return PyArray_Type.tp_as_number->nb_@oper@(a,b);
        case -2:
            /* use default handling */
            if (PyErr_Occurred()) {
                return NULL;
            }
            return PyGenericArrType_Type.tp_as_number->nb_@oper@(a,b);
        case -3:
            /*
             * special case for longdouble and clongdouble
             * because they have a recursive getitem in their dtype
             */
            Py_INCREF(Py_NotImplemented);
            return Py_NotImplemented;
    }

    /*
     * here we do the actual calculation with arg1 and arg2
     * as a function call.
     */
    PyErr_SetString(PyExc_TypeError, "complex long doubles do not support remainder");
    return NULL;
}

/**end repeat**/

/**begin repeat
 *
 * #name = half, float, double, longdouble, cfloat, cdouble, clongdouble#
 *
 */

/**begin repeat1
 *
 * #oper = lshift, rshift, and, or, xor#
 *
 */

#define @name@_@oper@ NULL

/**end repeat1**/

/**end repeat**/

/**begin repeat
 * #name = (byte, ubyte, short, ushort, int, uint,
 *             long, ulong, longlong, ulonglong,
 *             half, float, double, longdouble,
 *             cfloat, cdouble, clongdouble)*3,
 *
 *         byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong#
 *
 * #type = (npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *             npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *             npy_half, npy_float, npy_double, npy_longdouble,
 *             npy_cfloat, npy_cdouble, npy_clongdouble)*3,
 *
 *         npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong#
 *
 * #otype = (npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *             npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *             npy_half, npy_float, npy_double, npy_longdouble,
 *             npy_cfloat, npy_cdouble, npy_clongdouble)*2,
 *         npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *         npy_half, npy_float, npy_double, npy_longdouble,
 *         npy_float, npy_double, npy_longdouble,
 *
 *         npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 *         npy_long, npy_ulong, npy_longlong, npy_ulonglong#
 *
 * #OName = (Byte, UByte, Short, UShort, Int, UInt,
 *              Long, ULong, LongLong, ULongLong,
 *              Half, Float, Double, LongDouble,
 *              CFloat, CDouble, CLongDouble)*2,
 *          Byte, UByte, Short, UShort, Int, UInt,
 *          Long, ULong, LongLong, ULongLong,
 *          Half, Float, Double, LongDouble,
 *          Float, Double, LongDouble,
 *
 *          Byte, UByte, Short, UShort, Int, UInt,
 *          Long, ULong, LongLong, ULongLong#
 *
 * #oper = negative*17, positive*17, absolute*17, invert*10#
 */
static PyObject *
@name@_@oper@(PyObject *a)
{
    @type@ arg1;
    @otype@ out;
    PyObject *ret;

    switch(_@name@_convert_to_ctype(a, &arg1)) {
    case 0:
        break;
    case -1:
        /* can't cast both safely use different add function */
        Py_INCREF(Py_NotImplemented);
        return Py_NotImplemented;
    case -2:
        /* use default handling */
        if (PyErr_Occurred()) {
            return NULL;
        }
        return PyGenericArrType_Type.tp_as_number->nb_@oper@(a);
    }

    /*
     * here we do the actual calculation with arg1 and arg2
     * make it a function call.
     */

    @name@_ctype_@oper@(arg1, &out);

    ret = PyArrayScalar_New(@OName@);
    PyArrayScalar_ASSIGN(ret, @OName@, out);

    return ret;
}
/**end repeat**/

/**begin repeat
 *
 * #name = half, float, double, longdouble, cfloat, cdouble, clongdouble#
 */

#define @name@_invert NULL

/**end repeat**/

#define _IS_NONZERO(x) (x != 0)
/**begin repeat
 *
 * #name = byte, ubyte, short, ushort, int,
 *         uint, long, ulong, longlong, ulonglong,
 *         half, float, double, longdouble,
 *         cfloat, cdouble, clongdouble#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int,
 *         npy_uint, npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *         npy_half, npy_float, npy_double, npy_longdouble,
 *         npy_cfloat, npy_cdouble, npy_clongdouble#
 * #simp = 1*14, 0*3#
 * #nonzero = _IS_NONZERO*10, !npy_half_iszero, _IS_NONZERO*6#
 */
static int
@name@_bool(PyObject *a)
{
    int ret;
    @type@ arg1;

    if (_@name@_convert_to_ctype(a, &arg1) < 0) {
        if (PyErr_Occurred()) {
            return -1;
        }
        return PyGenericArrType_Type.tp_as_number->nb_bool(a);
    }

    /*
     * here we do the actual calculation with arg1 and arg2
     * make it a function call.
     */

#if @simp@
    ret = @nonzero@(arg1);
#else
    ret = (@nonzero@(arg1.real) || @nonzero@(arg1.imag));
#endif

    return ret;
}
/**end repeat**/
#undef _IS_NONZERO


static int
emit_complexwarning(void)
{
    static PyObject *cls = NULL;
    npy_cache_import("numpy.core", "ComplexWarning", &cls);
    if (cls == NULL) {
        return -1;
    }
    return PyErr_WarnEx(cls,
            "Casting complex values to real discards the imaginary part", 1);
}

/**begin repeat
 *
 * #name = byte, ubyte, short, ushort, int,
 *         uint, long, ulong, longlong, ulonglong,
 *         half, float, double, longdouble,
 *         cfloat, cdouble, clongdouble#
 *
 * #Name = Byte, UByte, Short, UShort, Int,
 *         UInt, Long, ULong, LongLong, ULongLong,
 *         Half, Float, Double, LongDouble,
 *         CFloat, CDouble, CLongDouble#
 *
 * #cmplx = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1#
 * #sign = (signed, unsigned)*5, , , , , , , #
 * #ctype = long*8, PY_LONG_LONG*2,
 *          double*3, npy_longdouble, double*2, npy_longdouble#
 * #to_ctype = , , , , , , , , , , npy_half_to_double, , , , , , #
 * #func = (PyLong_FromLong, PyLong_FromUnsignedLong)*4,
 *         PyLong_FromLongLong, PyLong_FromUnsignedLongLong,
 *         PyLong_FromDouble*3, npy_longdouble_to_PyLong,
 *         PyLong_FromDouble*2, npy_longdouble_to_PyLong#
 */
static PyObject *
@name@_int(PyObject *obj)
{
    PyObject *long_result;

#if @cmplx@
    @sign@ @ctype@ x = @to_ctype@(PyArrayScalar_VAL(obj, @Name@).real);
#else
    @sign@ @ctype@ x = @to_ctype@(PyArrayScalar_VAL(obj, @Name@));
#endif

#if @cmplx@
    if (emit_complexwarning() < 0) {
        return NULL;
    }
#endif

    long_result = @func@(x);
    if (long_result == NULL){
        return NULL;
    }

    return long_result;
}
/**end repeat**/

/**begin repeat
 *
 * #name = byte, ubyte, short, ushort, int, uint,
 *             long, ulong, longlong, ulonglong,
 *             half, float, double, longdouble,
 *             cfloat, cdouble, clongdouble#
 * #Name = Byte, UByte, Short, UShort, Int, UInt,
 *             Long, ULong, LongLong, ULongLong,
 *             Half, Float, Double, LongDouble,
 *             CFloat, CDouble, CLongDouble#
 * #cmplx = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1#
 * #to_ctype = , , , , , , , , , , npy_half_to_double, , , , , , #
 * #func = PyFloat_FromDouble*17#
 */
static NPY_INLINE PyObject *
@name@_float(PyObject *obj)
{
#if @cmplx@
    if (emit_complexwarning() < 0) {
        return NULL;
    }
    return @func@(@to_ctype@(PyArrayScalar_VAL(obj, @Name@).real));
#else
    return @func@(@to_ctype@(PyArrayScalar_VAL(obj, @Name@)));
#endif
}
/**end repeat**/

/**begin repeat
 * #oper = le, ge, lt, gt, eq, ne#
 * #op = <=, >=, <, >, ==, !=#
 * #halfop = npy_half_le, npy_half_ge, npy_half_lt,
 *           npy_half_gt, npy_half_eq, npy_half_ne#
 */
#define def_cmp_@oper@(arg1, arg2) (arg1 @op@ arg2)
#define cmplx_cmp_@oper@(arg1, arg2) ((arg1.real == arg2.real) ?        \
                                      arg1.imag @op@ arg2.imag :        \
                                      arg1.real @op@ arg2.real)
#define def_half_cmp_@oper@(arg1, arg2) @halfop@(arg1, arg2)
/**end repeat**/

/**begin repeat
 * #name = byte, ubyte, short, ushort, int, uint,
 *         long, ulong, longlong, ulonglong,
 *         half, float, double, longdouble,
 *         cfloat, cdouble, clongdouble#
 * #simp = def*10, def_half, def*3, cmplx*3#
 */
static PyObject*
@name@_richcompare(PyObject *self, PyObject *other, int cmp_op)
{
    npy_@name@ arg1, arg2;
    int out=0;

    RICHCMP_GIVE_UP_IF_NEEDED(self, other);

    switch(_@name@_convert2_to_ctypes(self, &arg1, other, &arg2)) {
    case 0:
        break;
    case -1:
        /* can't cast both safely use different add function */
    case -2:
        /* use ufunc */
        if (PyErr_Occurred()) {
            return NULL;
        }
        return PyGenericArrType_Type.tp_richcompare(self, other, cmp_op);
    case -3:
        /*
         * special case for longdouble and clongdouble
         * because they have a recursive getitem in their dtype
         */
        Py_INCREF(Py_NotImplemented);
        return Py_NotImplemented;
    }

    /* here we do the actual calculation with arg1 and arg2 */
    switch (cmp_op) {
    case Py_EQ:
        out = @simp@_cmp_eq(arg1, arg2);
        break;
    case Py_NE:
        out = @simp@_cmp_ne(arg1, arg2);
        break;
    case Py_LE:
        out = @simp@_cmp_le(arg1, arg2);
        break;
    case Py_GE:
        out = @simp@_cmp_ge(arg1, arg2);
        break;
    case Py_LT:
        out = @simp@_cmp_lt(arg1, arg2);
        break;
    case Py_GT:
        out = @simp@_cmp_gt(arg1, arg2);
        break;
    }

    if (out) {
        PyArrayScalar_RETURN_TRUE;
    }
    else {
        PyArrayScalar_RETURN_FALSE;
    }
}
/**end repeat**/

/**begin repeat
 *  #name = byte, ubyte, short, ushort, int, uint,
 *          long, ulong, longlong, ulonglong,
 *          half, float, double, longdouble,
 *          cfloat, cdouble, clongdouble#
**/
static PyNumberMethods @name@_as_number = {
    .nb_add = (binaryfunc)@name@_add,
    .nb_subtract = (binaryfunc)@name@_subtract,
    .nb_multiply = (binaryfunc)@name@_multiply,
    .nb_remainder = (binaryfunc)@name@_remainder,
    .nb_divmod = (binaryfunc)@name@_divmod,
    .nb_power = (ternaryfunc)@name@_power,
    .nb_negative = (unaryfunc)@name@_negative,
    .nb_positive = (unaryfunc)@name@_positive,
    .nb_absolute = (unaryfunc)@name@_absolute,
    .nb_bool = (inquiry)@name@_bool,
    .nb_invert = (unaryfunc)@name@_invert,
    .nb_lshift = (binaryfunc)@name@_lshift,
    .nb_rshift = (binaryfunc)@name@_rshift,
    .nb_and = (binaryfunc)@name@_and,
    .nb_xor = (binaryfunc)@name@_xor,
    .nb_or = (binaryfunc)@name@_or,
    .nb_int = (unaryfunc)@name@_int,
    .nb_float = (unaryfunc)@name@_float,
    .nb_floor_divide = (binaryfunc)@name@_floor_divide,
    .nb_true_divide = (binaryfunc)@name@_true_divide,
    /* TODO: This struct/initialization should not be split between files */
    .nb_index = (unaryfunc)NULL,  /* set in add_scalarmath below */
};
/**end repeat**/

NPY_NO_EXPORT void
add_scalarmath(void)
{
    /**begin repeat
     *  #name = byte, ubyte, short, ushort, int, uint,
     *          long, ulong, longlong, ulonglong,
     *          half, float, double, longdouble,
     *          cfloat, cdouble, clongdouble#
     *  #NAME = Byte, UByte, Short, UShort, Int, UInt,
     *          Long, ULong, LongLong, ULongLong,
     *          Half, Float, Double, LongDouble,
     *          CFloat, CDouble, CLongDouble#
     **/
    @name@_as_number.nb_index = Py@NAME@ArrType_Type.tp_as_number->nb_index;
    Py@NAME@ArrType_Type.tp_as_number = &(@name@_as_number);
    Py@NAME@ArrType_Type.tp_richcompare = @name@_richcompare;
    /**end repeat**/
}


NPY_NO_EXPORT int initscalarmath(PyObject * m)
{
    add_scalarmath();

    return 0;
}
